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ABSTRACT 


The use of computational fluid dynamics in scramjet engine component development is widespread in the 
existing literature. Unfortunately, the quantification of model-form uncertainties is rarely addressed with 
anything other than sensitivity studies, requiring that the computational results be intimately tied to and 
calibrated against existing test data. This practice must be replaced with a formal uncertainty quantification 
process for computational fluid dynamics to play an expanded role in the system design, development, and 
flight certification process. Due to ground test facility limitations, this expanded role is believed to be a 
requirement by some in the test and evaluation community if scramjet engines are to be given serious con- 
sideration as a viable propulsion device. An effort has been initiated at the NASA Langley Research Center 
to validate several turbulence closure models used for Reynolds-averaged simulations of scramjet isolator 
flows. The turbulence models considered were the Menter BSL, Menter SST, Wilcox 1998, Wilcox 2006, 
and the Gatski-Speziale explicit algebraic Reynolds stress models. The simulations were carried out using 
the VULCAN computational fluid dynamics package developed at the NASA Langley Research Center. A 
procedure to quantify the numerical errors was developed to account for discretization errors in the val- 
idation process. This procedure utilized the grid convergence index defined by Roache as a bounding 
estimate for the numerical error. The validation data was collected from a mechanically back-pressured 
constant area (1x2 inch) isolator model with an isolator entrance Mach number of 2.5. As expected, 
the model-form uncertainty was substantial for the shock-dominated, massively separated flowfield within 
the isolator as evidenced by a 6 duct height variation in shock train length depending on the turbulence 
model employed. Generally speaking, the turbulence models that did not include an explicit stress limiter 
more closely matched the measured surface pressures. This observation is somewhat surprising, given 
that stress-limiting models have generally been developed to better predict shock-separated flows. All of 
the models considered also failed to properly predict the shape and extent of the separated flow region 
caused by the shock boundary layer interactions. However, the best performing models were able to predict 
the isolator shock train length (an important metric for isolator operability margin) to within 1 isolator duct 
height. 

f Approved for public release; distribution is unlimited. 



INTRODUCTION 


A dual-mode scramjet engine (see Fig. 1 ) is composed of four main components: inlet, isolator, combustor, 
and exit nozzle. The inlet captures the air mass for the engine and conditions (compresses) the flow for sub- 
sequent combustion. The role of the isolator is to separate the combustion effects (i.e. combustion-induced 
pressure rise) from the inlet, providing sufficient margin to prevent engine un-start. The combustor contains 
the fuel injection devices through which fuel is introduced and subsequently mixed with the incoming air 
allowing for combustion. Finally, the nozzle expands the resulting exhaust gases to produce thrust for the 
vehicle. At low hypersonic Mach numbers, the heat release associated with the combustion process can 
lead to (or at least approach) thermal choking of the engine flow, resulting in elevated pressure levels in 
the combustor. The adverse pressure gradient that forms under these conditions leads to the formation of 
a pre-combustion shock train within the isolator. This flow is characterized by regions of massive flow sep- 
aration due to the adverse pressure gradient and the multiple shock/boundary layer interactions that result 
from this shock system. The isolator length required to contain the pressure rise induced by the combustion 
processes is of critical importance to the engine designer, since this drives the size (and weight) of the 
isolator section. 



Figure 1 : Dual-mode scramjet engine flowfield characteristics 

Computational Fluid Dynamics (CFD) has been used extensively to fill critical gaps present in scramjet 
engine research and development. The heavy reliance on CFD is driven by several factors that include: 


• difficulties with replicating hypersonic conditions in ground test facilities for sufficient periods of time 

• obstacles with utilizing advanced diagnostics in the harsh scramjet flow environments 

• lack of facilities large enough to test full-scale engines, other than possibly missile applications 


Reynolds-Averaged Simulations (RAS) have been used almost exclusively to fill these gaps due to the enor- 
mous computational resources required to perform scale-resolved calculations such as Large Eddy Simu- 
lation (LES). The epistemic uncertainties (i.e. systematic uncertainties) associated with imperfect physical 
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models used within a RAS framework are often found to dominate the overall uncertainty in simulations of 
scramjet flowpaths. 12 If CFD is to be used as a predictive tool in the design and analysis of scramjet com- 
ponents, these model-form uncertainties must be well quantified. The quantification of these uncertainties 
will require that the models in question be validated for their intended purposes (i.e. for realistic geometries 
in representative flow environments). To date, while the use of CFD is prevalent, very few efforts have been 
undertaken that truly attempt to validate models for component level simulations. Instead, the current state- 
of-the-art relies heavily on the experience of the CFD practitioner to estimate the model-form uncertainty 
associated with their simulations through simple sensitivity studies. This practice will have to be replaced 
with a formal uncertainty quantification process if CFD is to play an expanded role in the Design Develop- 
ment Research and Engineering community, Development Test and Evaluation community, and ultimately 
certification for flight. This is especially true with hypersonic air-breathing propulsion systems due to the 
environment, scale, and duration limitations of ground test facilities. The work described here represents 
a step toward enabling a process for quantifying model-form uncertainties for the isolator component of a 
dual-mode scramjet engine. In particular, the purpose of this effort is to document the initial stages of a for- 
mal model validation process underway at the NASA Langley Research Center for the flow physics present 
in scramjet isolators. 


FACILITY DESCRIPTION 


The Isolator Dynamics Research Laboratory (IDRL) has recently been brought on-line by the Flypersonic 
Airbreathing Propulsion Branch at the NASA Langley Research Center. A photograph and schematic of 
the facility hardware are given in Fig. 2. The flowpath is mounted vertically with facility flow directed from 
bottom to top. A unique aspect of this facility is the ability to translate and/or rotate the rig during facility 
operation. This capability was desired to simplify the application of advanced laser-based diagnostics for 
flowfield interrogation. The ability to move the test apparatus allows the diagnostic equipment to be set-up 
and aligned only once, so that the diagnostic components never have to be relocated. Air is supplied via 4 
flexible hoses to a settling chamber that houses a honeycomb flow straightener to remove swirl and reduce 
freestream turbulence levels. The settling chamber connects to a planar Mach 2.5 nozzle that was designed 
using an inviscid method-of-characteristics code coupled with a two-dimensional full Navier-Stokes solver. 
The 24-inch isolator test section has a height ( h ) of 1-inch and a width (w) of 2-inches, and is constructed 
through any combination of glass and metallic side walls. Flence, both wall measurements and laser-based 
in-stream measurements can be obtained simultaneously. The metal walls were installed for this effort, 
allowing for a total of 213 static pressure measurements in the arrangement shown in Fig. 3. The flow exits 
the isolator test section and enters a 15-inch long cylindrical expansion section that is machined with a 45° 
chamfer on the downstream end to mate with a back-pressure plug. The back-pressure plug is mounted in 
the spool section to a linear actuator that positions the plug with a precision of 1/1000 of an inch. The back- 
pressure plug simulates the thermal choking effect caused by combustion in a dual-mode scramjet engine, 
providing a means of creating the back pressured isolator flowfield of interest. A detailed discussion of the 
facility is given in Ref. 3. 


COMPUTATIONAL DESCRIPTION 


The computational domain considered spans the region from the facility nozzle plenum to the end of the 
24-inch isolator duct. The geometrical symmetry present in this flowpath, and the fact that only steady-state 
simulations are to be performed, permitted the consideration of only one quarter of the facility flowpath 
as shown in Fig. 4. A structured grid with a total of 85,196,800 cells was created for the simulations 
with 18,636,800 cells placed in the nozzle block (560 x 160 x 208), and the remaining 66,560,000 cells 
(2000 x 1 60 x 208) in the isolator duct. Two coarser grids (1 0,649,600 and 1 ,331 ,200 cells) were obtained 
from this parent grid via a factor of 2 coarsening in each coordinate direction to allow for a formal assessment 
of the numerical errors associated with the simulations. The grid points were clustered near the walls such 
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Figure 2: Photograph (left) and schematic (right) of the Isolator Dynamics Research Laboratory hardware 
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Figure 3: Instrumentation layout for the metallic isolator side walls: static pressure ports (blue), thermocou- 
ples (red), high frequency pressure ports (green) 


that the maximum y + value at any cell center adjacent to the surface was no greater than 1 (approximately) 
for the coarsest grid considered. The grid was also clustered in the streamwise direction at the nozzle throat 
(Ax = 0.0025 h) with the spacing gradually increased to 0.01 h at the nozzle exit. The streamwise spacing 
was held fixed at 0.01 h for the front half of the isolator to resolve the region where the shock train was 
expected to reside, and was stretched in the second half of the isolator to a value of 0.025 h at the isolator 
exit. The maximum grid spacing in the vertical (y) and spanwise (z) directions was 0.01 h (adjacent to the 
planes of symmetry). The entire grid was split up into a total of 1280 grid blocks which yielded good (95% 
or better) load balance statistics for up to 438 processors. 

All computational results were obtained using the VULCAN (Viscous Upwind algorithm for Complex flow 
ANalysis) Navier-Stokes code. 4 This CFD package is developed and maintained by researchers in the 
Hypersonic Airbreathing Propulsion Branch at the NASA Langley Research Center, and the development 
of this software has predominantly been driven by the needs of the scramjet community. The CFD data 
obtained in this effort were acquired by integrating the Reynolds-Averaged Navier-Stokes equations until 
steady-state conditions were achieved. All of the solutions were advanced in pseudo-time via an incom- 
plete LU factorization scheme (with planar relaxation) using a Courant-Friedrichs-Lewy (CFL) number of 
100. The inviscid fluxes were evaluated using the Low-Diffusion Flux Splitting Scheme of Edwards, 5 with 
cell interface variable reconstruction achieved via the k= 1/3 Monotone Upstream-centered Scheme for Con- 
servation Laws. The Koren flux limiter 6 was utilized to avoid spurious oscillations during this reconstruction 
process. The viscous fluxes were evaluated using 2 nd -or6e\ accurate central differences with the viscosity of 
air computed from the polynomial fit of McBride. 7,8 The molecular and turbulent Prandtl numbers were set 
to 0.72 and 0.9, respectively. An adiabatic no-slip condition was applied at all solid surfaces. The subsonic 
nozzle plenum inflow condition was specified based on the measured plenum pressure and temperature, 
with Mach number extrapolated from the interior. Finally, all properties were extrapolated at the subsonic 
exit of the isolator with the exception of the static pressure. The static pressure was set to the value ob- 
tained by averaging the 20 surface pressure measurements taken around the periphery of the isolator duct 
at the x=22 h station. This station was highly instrumented with pressure taps for the specific purpose of 
providing a well characterized outflow boundary condition. 
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Figure 4: Grid created for the nozzle/isolator assembly (every other point removed for clarity): nozzle section 
(top), first half of isolator section (middle), second half of isolator section (bottom) 
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Five different turbulence models were assessed in this effort: Menter BSL, 9 Menter SST, 9 Wilcox 1998, 10 
Wilcox 2006, 11 and the Gatski-Speziale Explicit Algebraic Stress Model (EASM). 12 The implementation of 
these models is precisely as described in the references cited above, with the following (potential) excep- 
tions: 


• The “exact” production term in the transport equation for the turbulence kinetic energy (i.e. t u ) was 
consistently utilized. 

• The production of turbulence kinetic energy (k) was limited to be no more than 1 0 times the destruction 
of k to discourage spurious production of turbulence across shock waves. 

• The surface boundary condition for the specific dissipation rate ( co ) was set as suggested by Menter 9 
to be 10 times the derived asymptotic value. 

• The modification to the destruction term coefficient for the co equation in the Wilcox models that ad- 
dresses the “round jet / plane jet anomaly” was omitted due to undesirable (and unintended) effects 
noted in the shock separated regions. 


It is also worth noting that a variety of sensitivity studies were performed to address the impact of modeling 
decisions other than the choice of turbulence model. These simulations were performed to ensure that any 
variability in the simulation results of interest were small relative to variations expected as the turbulence 
models are varied. Modeling sensitivities that were considered include: 


• choice of turbulent Prandtl number (varied between 0.7 and 1.1) 

• choice of model for the turbulence production term (exact or approximate formulation based on the 
vorticity magnitude) 

• choice of limiting value placed on the production to destruction ratio for k (limited by 10 or unlimited) 

• choice of viscous constitutive relation (Sutherland law or McBride polynomial) 


Although not detailed in this document, the solution variability associated with these sensitivities were small 
relative to variations noted when changing the turbulence model; confirming that the turbulence model 
closure is the dominant source of epistemic uncertainty for this flowfield class. Finally, additional simulations 
were performed to check for hysteresis since it is known that flows with massive separation can be sensitive 
to how the separations develop during the process of converging to a steady-state. In particular, simulations 
were performed using an initial condition based on an existing converged solution (obtained from the use of 
a different turbulence model) to ensure that the results matched those obtained by the default initialization 
procedure. While the authors do not definitively claim that the solutions obtained are not in some way 
influenced by how the flowfields were initialized, these checks at least ensured that the solutions obtained 
were achievable via more than one numerical solution path. 


RESULTS AND DISCUSSION 


The IDRL facility was designed from its inception with CFD validation in mind. Toward this end, both 
Experimental Fluid Dynamics (EFD) and CFD researchers were involved in all stages of the design process. 
This collaboration was crucial to ensure that all data requirements (as defined by the CFD team) and 
measurement capabilities (as defined by the EFD team) were communicated. This communication led to 
the determination of the diagnostic approaches to employ, and the hierarchy of measurements required, to 
maximize the value of the experimental data. Additional highlights from this collaboration are listed below 
(the details of which can be found in Ref. 13): 
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• CFD simulations were conducted to define facility requirements, including: 

- sensitivity studies to determine the control authority (precision) required for the back-pressure 
plug 

- sensitivity studies to determine the control authority requirements for the facility supply pressure 

- sensitivity studies to determine the influence of supply temperature (and surface temperature) to 
determine whether a temperature control system should be considered 

- sensitivity studies to determine the influence of facility supply air turbulence levels 

• Instrumentation layout was defined (to a large extent) by the CFD team with particular emphasis 
placed on measurements required to specify and/or assess boundary conditions. 

• CFD simulations were performed to determine the ramification of slot openings present in the specially 
designed isolator side walls required by some of the advanced laser-based diagnostics (not present 
in this study). 


Solution Verification 


The initial step of any validation effort requires an assessment of numerical errors. Numerical errors are 
denoted here as spatial discretization errors and iterative convergence errors. Temporal discretization er- 
rors are of no consequence for this study since only steady-state simulations are considered. The most 
compelling proof of iterative convergence is the situation where the residual norms have been reduced to 
machine accuracy. Unfortunately, this situation is rarely realized in practice, particularly when simulating 
complex flows with shock waves and regions of massive flow separation. Alternatively, the L 2 norm of the 
residual is monitored until it levels out about some acceptably low value. In addition, the variability of integral 
properties of interest such as mass flow rates and surfaces loads are checked to offer further assurance that 
an adequate level of convergence has been achieved. In general, the following iterative convergence state- 
ments were satisfied for each simulation discussed in this effort (also see Fig. 5 for a typical convergence 
history): 


• The value of the residual L 2 norm was reduced by at least 5 orders of magnitude. 

• The integrated surface load time histories were unchanged (to at least 5 significant digits) over the 
last 5000 iteration cycles. 

• The integrated mass flow rate was constant at every streamwise grid plane to 4 significant digits. 

The Grid Convergence Index (GCI) 14 was used to quantify the spatial discretization (i.e. grid convergence) 
error. The GCI is a grid convergence estimator derived from the generalized Richardson Extrapolation 
formula, and can be written as follows: 

GCI = F s — (1) 
rP- 1 

In this expression, e is the difference of some functional “/” evaluated using two different grid resolutions 
with refinement ratio r, i.e. 

e = fi-.n ( 2 ) 

and p is either the observed order of accuracy of the numerical scheme or taken to be 2 (i.e. 2" rf -order). 
Finally, F s is a “factor of safety” with recommended values taken to be either 3 (if the observed order 
of accuracy is assumed to be the theoretical value) or 1 .25 (if the observed order of accuracy has been 
rigorously determined). It should be emphasized that the GCI, while based on Richardson Extrapolation, 
is not meant to be a “best estimate” of the numerical error. Instead, the intent is to provide a reasonable 
bound on the discretization error. 
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Figure 5: Iterative convergence measures (circles show start of coarse, medium, and fine grid sequences) 


The observed order of accuracy was extracted from the isolator simulations performed in this effort. A 
minimum of three grids are required to facilitate the extraction of the observed order of accuracy via the 
following expression: 

p ’' n [ f 0 ) ' (3) 

The coarse grids were obtained by systematically removing every other point from the grid that was gen- 
erated, resulting in a refinement (or coarsening) ratio of r= 2. The subscript of each functional denotes the 
grid level with “1” representing the finest grid and “3” representing the coarsest. Note that this expression 
implicitly assumes that solutions obtained from each grid level are within the asymptotic range of conver- 
gence. There is no way to formally prove that each grid level satisfies this constraint (without adding yet 
another grid level). However, the coarsest grid was designed to maintain a y + value near unity (or less) 
adjacent to each solid surface, and was judged to offer adequate resolution based on past experience with 
performing simulations of this type. Moreover, sanity checks were performed using a variety of functionals 
to check that the order of accuracy extracted was consistent between functionals of a given class. 

The conditions chosen to assess the observed order of accuracy were based on flow conditions that result 
in a pressure ratio across the isolator shock train that is approximately 80% of the normal shock value based 
on the approach flow Mach number. This condition represents a realistic operability limit that dual-mode 
scramjets are often designed to, and is a close match to the conditions that will be utilized for the validation 
exercises. The precise conditions chosen are given in Table 1 . The turbulence model selected for this 
exercise was the Menter BSL model. 

The isolator flowfield obtained from these flow conditions is presented in Fig. 6. The global isolator shock 
train structure (visualized by the velocity divergence) consists of leading oblique shock waves followed by 
alternating expansions and quasi-normal shocks. The top image shows static pressure contours super- 
imposed on an iso-surface of the total pressure (P a = 300 kPa). This image provides a visual of the core flow 
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Table 1 : Facility flow conditions chosen for order of accuracy assessment 


Nozzle Plenum Pressure [kPa] 

861.85 

Nozzle Plenum Temperature [K] 

298.15 

Nozzle Exit Mach Number 

2.5* 

Nozzle Exit Pressure [kPa] 

50.44 * 

Isolator Exit Pressure [kPa] 

296.50 


* Nominal value computed from isentropic relationships 


displacement caused by the interactions between the shock waves and the turbulent boundary layer that 
engulfs the isolator walls. The area blockage caused by the separated flow in the corners is a dominant 
feature that dictates (to a large extent) the shock train flow structure. The Mach number contours show 
a boundary layer structure on the vertical side walls that is notably thicker near the center of the side 
wall. This feature is a result of the expansion processes in the facility nozzle, and makes this particular 
boundary layer more susceptible to large-scale flow separation. Mach disks are also seen in the center 
of the duct, and while it is difficult to visualize in this image, the interaction of the leading oblique shock 
wave and the Mach disk results in the formation of a shear layer. The flow in the isolator initially separates 
in the corners of the duct where the shear forces are weakest. This separation results in a conical shock 
that spreads diagonally across the isolator; diverting the flow toward the center of the duct. The high 
pressure that is produced in the corners also drives the the boundary layer fluid away from the corners 
along the adjacent walls. The side walls (z=constant vertical surfaces) are half the width of the top and 
bottom walls (y=constant spanwise surfaces), so the effect of the adverse pressure gradient formed in the 
corners influences the center of the vertical walls more quickly than their spanwise counterparts. The 
additional width offered by the spanwise surfaces allows sufficient space for a second separation region to 
form at the z=0 centerplane of the duct. This separation bubble forms a strong oblique shock followed by 
a strong expansion as this separation zone re-attaches to the surface. The resulting flow remains attached 
along the isolator centerplane even after further shock reflections downstream. The separation zones that 
originate in the corners, however, gradually extend across the side walls and persist for several duct heights 
further downstream. Moreover, the increased turbulence levels present in these separated flow zones do 
not dissipate as the flow progresses downstream. The eventual result of this unbalanced production of 
turbulence is a very large “effective” viscosity across the entire width of the duct that is 20,000 times larger 
than the local fluid viscosity. The net effect is a flowfield that exits the isolator with properties representative 
of a Reynolds number that is 4 orders of magnitude lower than the flow entering the isolator. This behavior 
is the mechanism that allows a Reynolds-averaged simulation to produce a steady-state solution for this 
flow. 

A total of 5 functionals (representing 2 functional classes) were extracted from this flowfield to perform the 
order of accuracy assessment. The surface shear (skin friction) at the nozzle exit defined one functional 
class of interest. This parameter provides a measure of how resistant the approach flow boundary layer is 
to flow separation. The wall shear along the center of both the vertical and spanwise surfaces was extracted 
and used to compute the order of accuracy for this functional class. The second functional class considered 
involved measures of isolator shock strength. Local measures of shock properties were avoided since the 
theory behind Richardson Extrapolation (which is based on a Taylor series expansion) breaks down for 
properties that are discontinuous. Instead, the integrated isolator surface force components and incipient 
shock locations were considered for the order of accuracy assessment within this class. 

The results of the order of accuracy assessment are given in Table 2, showing that the observed order of 
accuracy is only 1 "-order even though a nominally 2 nd -or6er spatial scheme was utilized. This behavior, 
while troubling, can be explained by noting that a non-linear TVD (Total Variation Diminishing) flux limiter 
was invoked to capture the anticipated strong shocks with minimal dispersion. The TVD property of these 
limiters forces the numerical scheme to become locally 1 "-order accurate whenever a flow property extrema 
is encountered. As pointed out by Roache, 14 it is not uncommon for errors introduced by localized 1 "-order 
treatments to propagate to other areas of the flowfield. This is particularly true for systems of equations 
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Figure 6: Isolator flowfield characteristics: (a) total pressure iso-surface colored by static pressure con- 
tours, (b) Mach number contours (c) velocity divergence contours, (d) turbulent to molecular viscosity ratio 
contours 


where more than one family of characteristics are present for information propagation. It should also be 
emphasized that although the order of accuracy assessment shows the rate of grid convergence to be 1 st - 
order, the solution accuracy is considerably higher than what would have been obtained using a globally 
I s ' -order algorithm on the same grid. In other words, the extracted order of accuracy simply provides a 
statement on the rate of grid convergence, and in isolation should not be used to judge solution accuracy. 


Table 2: Extracted order of accuracy for each functional 


Functional 

Fine Grid 

Medium Grid 

Coarse Grid 

Accuracy Order 

Nozzle Exit C f (spanwise centerline) 

1 .0956xl0 -4 

1.0831xl0 -4 

1.0572xl0- 4 

1.05 

Nozzle Exit Cf (vertical centerline) 

1 .0042xl0 -4 

9.9776xl0 -5 

9.8142xl0 -5 

1.33 

Isolator Surface F x [N] 

6.1311 

6.0661 

5.9228 

1.14 

Isolator Surface F y [N] 

2735.4 

2698.2 

2625.3 

0.97 

Isolator Surface F z [N] 

1368.2 

1349.6 

1312.9 

0.98 

Spanwise centerline shock location * 

8.8734 

9.1631 

9.7621 

1.05 

Vertical centerline shock location * 

8.2825 

8.5850 

9.2228 

1.08 

Corner shock location * 

7.7942 

8.0967 

8.7417 

1.09 


* Shock location defined as the distance from the nozzle exit (x/h) where the shock-induced pressure rise is 2.5 psi (17.237 kPa) 


Having established the order of accuracy of the simulations, the following recipe was utilized to quantify (in 
the form of an error band) the numerical error associated with the simulations for all subsequent validation 
exercises. First, it was assumed that the order of accuracy extracted from the simulations described above is 
applicable to the simulations performed for each of the turbulence models being validated. This was deemed 
to be a reasonable assumption given the fact that the turbulence models considered in this effort (with the 
exception of the EASM formulation) are functionally quite similar. Moreover, given that the observed order of 
accuracy was only 1 -order, the use of p = 1 in the expression for the GCI (see Eq. 1 ) provides a conservative 
estimate of the numerical error. The order of accuracy would only increase (if it changed at all) if this process 
were to be repeated for each turbulence model, so in the spirit of attempting to bound the error (rather than 
get a best estimate of it), the value of p was set to 1.0 for all GCI evaluations. This approach allows the 
use of only the medium and coarse grids in all remaining simulations for validation purposes; reducing 
the computational costs by at least an order of magnitude. Given the observed order of accuracy and the 
relationship for the GCI, an error band for the shock train position is computed for inclusion in the plots as 
an error bar. The error bar is one-sided, rather than symmetrically placed about the computed data point, 
because as the grid is resolved the variation of a given functional is monotonic (provided that each grid is 
within the asymptotic convergence range). 

An illustration of this strategy for reporting the numerical error is given in Fig. 7. In this example, the 
numerical error “bound” for the spanwise symmetry plane surface shock position has been quantified using 
the GCI for this functional given the observed order of accuracy {p= 1) and the coarse and medium grid 
shock position values taken from Table 2. The spanwise centerline surface pressure trace extracted from 
the fine grid solution has been included as a surrogate “truth” value. The error bars based on the GCI 
(computed for both the medium grid and coarse grid solutions) clearly bracket every data point of the fine 
grid simulation. While this does not prove that the numerical error has been rigorously bounded, it does 
provide some confidence that the approach is yielding a reasonable quantified bound of the numerical 
error. As a final note, Fig. 8 compares all three pressure traces, with the coarse and medium grid results 
translated to provide a best-fit to the fine grid data. The close match after translation supports the choice 
of using shock position as an error indicator rather than the discrete surface pressure values which vary 
discontinuously across shocks. The fact that the numerical error primarily resulted in a delayed initiation of 
the shock train, and did not alter the detailed shock structure to any meaningful degree, allowed the use 
of a scalar GCI value (shock position) rather than a more complicated approach based on an array of GCI 
values computed from spatially varying properties. 


12 






Figure 7: Comparison of fine grid wall pressure (spanwise centerline) with the medium grid results (left) and 
the coarse grid results (right); error bars represent the GCI estimated error bound 




Figure 8: Comparison of fine grid wall pressure (spanwise centerline) with the translated medium grid 
results (left) and the translated coarse grid results (right) 
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Model Validation 


The conditions of interest for the validation exercises are isolator conditions that result in a pressure ratio 
across the isolator shock train that approaches the operability bound of typical dual-mode scramjets. The 
IDRL isolator shock train length (L lS0 ) plotted as a percentage of the normal shock pressure ratio (based on 
the nominal isolator entrance Mach number of 2.5) is given in Fig. 9. The shock train length is defined as 
the distance from the location of the incipient shock-induced surface pressure rise to the exit of the isolator. 
This plot shows that the ability of an isolator to hold a given pressure ratio diminishes rapidly as the normal 
shock pressure ratio is approached. For the isolator conditions considered here, a 9% increase in the back 
pressure from 82.25% to 89.82% of the normal shock pressure rise changes L iso by 8.6 x/h (or 35% of the 
available isolator length). The sensitivity of the shock position to the back pressure value in this regime is 
clearly a challenge for any simulation effort to accurately predict. This sensitivity also explains why dual- 
mode scramjets have historically been designed to limit the total combustor-induced pressure rise during 
low-speed operation to approximately 80% of the normal shock value based on the isolator approach flow 
Mach number. 



Figure 9: Isolator shock train length as a function of the normal shock pressure rise for a Mach 2.5 flow 
(dashed portion of the line represents an operability range beyond that of current designs) 

The facility isolator conditions (with and without back pressure) considered for this validation exercise are 
given in Table 3. The measured conditions in this table, and the measured surface pressures in the plots 
that follow, are an ensemble average of measurements taken over an approximate 5 second test interval. 
The values in parenthesis represent the rms deviation in the measured values relative to the mean value 
obtained over the 5 second averaging interval. The IDRL facility has only very recently been brought on- 
line, and the uncertainty quantification process for the measurements has not been completed. As a result, 
only the measurement deviation rms over the test window is being reported at this time. 

A comparison of surface pressure obtained from the medium grid with measurements for the fully super- 
sonic (no applied back pressure) conditions is shown in Fig. 10. The simulated results shown are only those 
obtained with the Menter BSL turbulence model, since all of the turbulence models considered in this effort 
predict essentially the same surface pressure distribution when a back pressure is not applied. Two distinct 
observations are immediately apparent when comparing the simulations to the measured values. The first 
observation is that the residual wave structure exiting the nozzle is significantly stronger than what has 
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Table 3: Facility flow conditions selected for model validation 


Facility Condition 

Fully Supersonic 

Back Pressured 

Nozzle Plenum Pressure [kPa] 

860.30 (0.013) 

845.29 (0.025) 

Nozzle Plenum Temperature [K] 

297.00 (0.018) 

292.91 (0.023) 

Nozzle Exit Mach Number 

2.5* 

2.5* 

Nozzle Exit Pressure [kPa] 

50.35 * 

49.47 * 

Isolator Exit Pressure [kPa] 

N/A 

304.51 (0.097) 


* Nominal value computed from isentropic relationships 


been predicted. This result was expected to some extent, since there are known discrepancies between the 
as-designed (used for the simulations) versus the as-built nozzle geometry. The discrepancies are within 
the fabrication tolerances of the nozzle, but the team has not yet had the opportunity to address the impact 
of these geometrical differences. This will be a focal point in the next stage of this research effort. The 
more troubling observation is the systematic under-prediction of pressure throughout the entire length of 
the isolator. The two most obvious explanations for this particular disagreement are either a mismatch in 
the nozzle exit to throat area ratio or some sort of instrumentation problem. A careful examination of the 
nozzle assembly has ruled out the nozzle area ratio as a potential culprit. However, a comparison of the 
measured mass flow rate through the facility (via a Venturi flow meter) with the value obtained from choked 
flow relationships (given the facility plenum conditions and the nozzle throat area) revealed a discrepancy of 
approximately 3.3%. The measured mass flow rate was higher than the value obtained from the choked flow 
equations with a unity discharge coefficient, which clearly points to some issue with the instrumentation. 
Upon further examination, a problem with the facility plenum pressure measurement was uncovered, and 
further facility shake-down activities are underway to determine if the resolution of this problem eliminates 
the discrepancy with the Venturi flow rate measurement. 



Figure 10: Comparison of measured wall pressure with computed results (medium grid) obtained for condi- 
tions without an imposed back pressure 

Unfortunately the 3.3% discrepancy between measured and deduced flow rate, while in the right direction 
relative to the simulations, accounts for only about half of the difference between the simulated and mea- 
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sured wall pressure values. Fig. 10 shows something closer to a 7.5% increase in nozzle exit pressure is 
required to match the measured pressure levels in the isolator. As noted in Fig. 9, the shock train location is 
highly sensitive to the precise pressure rise across the shock train. Since the simulations use the measured 
back pressure as an outflow boundary condition, it is critical that the static pressure entering the isolator 
also be matched. Given that an issue was uncovered with the plenum pressure measurement, which is a 
measured value used as an inflow boundary condition, the decision was made to increase the plenum pres- 
sures given in Table 3 by 7.5% for the simulations. There are other boundary condition modifications that 
could have been made to match the overall pressure rise through the isolator, so an additional simulation 
was performed to assess this sensitivity. The additional simulation involved reducing the back pressure by 
7.5%, while keeping the nozzle plenum pressure at the measured value. The difference in the normalized 
wall pressure distributions from each simulation is illustrated in Fig. 1 1 . This difference, due solely to the 
slight disparity in Reynolds number, is smaller than the quantified numerical error, suggesting that the de- 
tailed method used to match the simulated conditions in to and out of the isolator is of little consequence. 
However, until the root cause of the discrepancy between measured and computed nozzle exit conditions is 
understood, the validation exercise will be incomplete. Hence, comparisons with the measured pressures 
that follow in the remainder of this paper should be regarded as preliminary results. 



Figure 1 1 : Normalized wall pressure distributions obtained via alternative boundary conditions that match 
the measured pressure rise 

A comparison of the back pressured results obtained using the modified plenum pressure for each turbu- 
lence model are shown in Fig. 12. The computed results shown are those obtained from simulations that 
used the medium grid. The images on the left compare the pressure along the z=0 centerline of the top and 
bottom walls, while the images on the right compare pressure values along the y=0 centerline of the side 
walls. The error bars associated with the CFD results represent the numerical discretization error bound 
for shock position based on the approach described in the solution verification section. Overall, there is at 
least a 6 duct height variation in shock position depending on the turbulence model employed, so clearly 
the model-form uncertainty is significant. Some general observations drawn from these comparisons are 
listed below: 


• The measured pressure distribution along the spanwise centerline (z=0) shows a monotonic rise, while 
the simulations do not. 
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• The computed side wall y=0 centerline pressure distributions display a monotonic rise (or at least 
nearly monotonic) which is consistent with the measurements. 

• The two turbulence models that have an explicit stress limiter (Menter SST and Wilcox 2006) are the 
worst performers for this shock-dominated massively separated flow. The stress limited models also 
show opposite trends with respect to shock train position. 

• To within the discretization error uncertainty, the remaining turbulence models predict a shock train 
length that is within 1 duct height of the measured value. 


These observations are discussed in further detail in the paragraphs that follow. 

Perhaps the most noticeable difference between the predicted and measured surface pressure distributions 
is the fact that the measured pressure distribution along the z= 0 centerline shows a monotonic rise, while 
the computed distribution does not. A popular explanation for this behavior 2 is the observation that highly 
back-pressured isolator flowfields often exhibit large scale unsteadiness with shock train motions that can 
alter the L iso at any given instance in time by several duct heights. Hence, the measured static surface 
pressure through the shock train actually represents some effective pressure gradient that results from the 
unsteady movement of the shock system. This is contrasted with a steady-state simulation where the shock 
system is fixed in time with no smoothing of the pressure rise due to unsteady shock motions. However, one 
can always ensemble average (in time) an unsteady system with shock waves to obtain a time-averaged 
state, and the resulting state will contain a stationary shock system. Hence, a less literal interpretation of 
this phenomenon revolves around what can reasonably be expected from a single point turbulence closure 
approximation. The closure models developed for Reynolds-averaged simulations are almost exclusively 
derived for attached shear-dominated flows. The ability of these models to accurately model the detailed 
effects caused by multiple unsteady shock waves interacting with an attached approach boundary layer is 
certainly questionable. In other words, the reason that the pressure rise is not monotonic may simply be 
a result of the fact that the shape and size of the recirculation zones are not accurately predicted by the 
models. Examples of efforts that attempt to model the complex effects of shock unsteadiness on Reynolds- 
averaged closure models are described in Refs. 15 and 16. 

To illustrate the impact of the separated flow structure on the surface pressure, a simulation was performed 
using a modified IDRL isolator geometry where the width was increased from 2-inches to 4-inches (resulting 
in a cross-sectional aspect ratio of 4 instead of 2). The reasoning behind considering this geometry is based 
on this author’s observation that steady-state Reynolds-averaged simulations tend to exaggerate the three- 
dimensional surface pressure features in the shock separated flow portion of the isolator. By increasing 
the aspect ratio of the isolator, the three-dimensional corner and side-wall effects are expected to be less 
dominant. A comparison of the surface pressure along the z=0 centerline for both isolator geometries is 
given in Fig. 13. Both of these simulations utilized the Menter BSL model. The increased width has resulted 
in a longer L iso , as well as a nearly monotonic rise in wall pressure. Hence, three-dimensional effects from 
the corners and side walls clearly have a significant impact on the isolator flowfield structure. The flow 
structure for both the 1 x 2 and 1 x 4 isolators are shown in Figs. 14 and 15, respectively. While there 
are several notable differences in the flows, the primary reason for the more monotonic centerline pressure 
rise for the larger aspect ratio case is the fact that the shock separated zone did not immediately re-attach 
along the z=0 centerline. The additional viscous damping made available by the larger separation region 
enhances diffusion; allowing the discrete pressure jumps across each shock reflection to smooth out prior 
to reaching the isolator surface. 

Another item that deserves attention is the peculiar observation noted with the Menter SST and Wilcox 2006 
models. Both of these models have built-in stress limiters that are designed to be more prone to separate 
when in the presence of an adverse pressure gradient. The Menter SST surface pressure values showed 
the expected trend of earlier boundary layer separation relative to the BSL version of the model. In fact, the 
SST model predicted the shock train to initiate 2 to 3 duct heights further upstream than the measurements 
indicate. Menter calibrated the SST model for transonic flows, and experience has shown 11 that it tends 
to over-predict the extent of separation for higher Mach number conditions. Hence, the over-prediction of 
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Pressure [kPa] Pressure [kPa] Pressure [kPa] 







Figure 12: Comparison of measured wall pressure with computed results (medium grid) obtained for each 
turbulence model: top and bottom walls (left), side walls (right) 
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Figure 13: Comparison of spanwise centerline (z= 0) surface pressure for the (1 x 2) and (1 x 4) geometries 

L is0 was to some extent expected for this application. The Wilcox 2006 model, on the other hand, appears 
to have discouraged flow separation as compared with the other models considered (even those without a 
stress limiter). This result was not anticipated. Upon comparing this model to the Menter SST formulation, 
one finds that there are two fundamental functional differences other than the values specified for the 
constants. The first difference is the scalar parameter that controls the activation of the stress limiter. The 
Menter SST formulation uses the magnitude of the vorticity vector, while the Wilcox 2006 formulation uses 
the magnitude of the strain rate tensor. The second difference involves the turbulent diffusion terms in the 
transport equations for the turbulence kinetic energy and specific dissipation rate. The Wilcox 2006 model 
uses the non-limited eddy viscosity for these terms, while the Menter SST model uses the stress-limited 
eddy viscosity consistently for both the Navier-Stokes and turbulence transport equations. 

A variant of the Wilcox 2006 model that uses the vorticity vector magnitude in lieu of the strain rate tensor 
magnitude is described in Ref. 1 1 . This version of the model is touted as an acceptable substitute for the 
strain rate version by the model developer. However, it should be noted that practically all of the test cases 
used for the development of the Wilcox 2006 model are two-dimensional problems where there are far less 
terms involved that differentiate the strain rate tensor magnitude from the vorticity vector magnitude. An 
additional simulation was performed with this variant of the Wilcox 2006 model to examine the impact that 
this modification has on the stress limiter behavior. A comparison of the z=0 centerline pressure distributions 
obtained from each variant of the Wilcox 2006 model is shown in Fig. 16. The Wilcox 2006 model that uses 
the vorticity-based stress limiter shows flow separation at an earlier streamwise station. However, the shock 
train position is still further downstream than what was predicted by Wilcox 1 998 model, which does not have 
a stress limiter. If one also alters the diffusion terms in the turbulent transport equations for k and a> to utilize 
the stress-limited eddy viscosity, one finds that this modified version of the Wilcox 2006 model predicts the 
same trend as the Menter SST model (see Fig. 16); in the sense that it predicts earlier separation than a 
similar model that does not invoke a stress limiter. The predicted shock system is further downstream than 
the Menter SST results because the Wilcox 2006 model uses a smaller stress limiter coefficient than the 
Menter SST model (to reduce the tendency to over-predict the extent of separation that plagues the SST 
model at higher Mach numbers). It should be emphasized that the intent of this exercise was to simply 
explain the differences in behavior between the Menter SST and Wilcox 2006 model and should not be 
interpreted as an endorsement for changing either of the models. 
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Figure 14: Isolator (1 x 2) flowfield characteristics: (a) total pressure iso-surface colored by static pressure 
contours, (b) Mach number contours (c) velocity divergence contours, (d) turbulent to molecular viscosity 
ratio contours 
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Figure 15: Isolator (1 x 4) flowfield characteristics: (a) total pressure iso-surface colored by static pressure 
contours, (b) Mach number contours (c) velocity divergence contours, (d) turbulent to molecular viscosity 
ratio contours 
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Figure 16: Effect of Wilcox 2006 modifications on the wall pressure (spanwise centerline); vorticity-based 
stress limiter (left), vorticity-based stress limiter and stress-limited diffusion in the turbulence transport equa- 
tions (right) 


SUMMARY AND FUTURE WORK 


A research effort has been initiated to validate several turbulence closure models used for Reynolds- 
averaged simulations of scramjet isolator flows. This effort was crafted from its inception to couple both 
experimental and computational research teams to maximize the effectiveness of the validation exercise. 
The turbulence models considered in this validation effort were the Menter BSL, Menter SST, Wilcox 1998, 
Wilcox 2006, and the Gatski-Speziale explicit algebraic Reynolds stress models. A procedure to quantify 
the numerical errors was developed to account for discretization errors in the validation process. This proce- 
dure utilized the grid convergence index defined by Roache as a bounding estimate for the numerical error. 
The validation data, which presently is limited to surface pressure measurements, was collected from a me- 
chanically back-pressured constant area (1x2 inch) isolator model with an isolator entrance Mach number 
of 2.5. Of the turbulence models considered, the models that did not include an explicit stress limiter more 
closely matched the measured surface pressures. This observation is somewhat surprising in a general 
sense, since stress limiters have generally been developed to better predict separation for flows subjected 
to an adverse pressure gradient. All of the models considered also failed to properly predict the shape 
and extent of the separated flow region caused by the shock boundary layer interactions. This conclusion 
was drawn via indirect evidence by noting that the measured surface pressures increased monotonically 
through the shock train, while the CFD predictions did not. Finally, while each of the turbulence models 
failed to properly predict all of the flow structure details associated with this massively separated flow, the 
best performing models were able to predict the shock train length to within 1 isolator duct height. The 
prediction of this metric is critically important, since it drives the overall length (and weight) of the isolator 
section. 

The validation of closure models is a process, and this effort has documented only the initial stages of this 
process. Future work will focus on quantifying the uncertainty of the experimental effort. This will not only 
include the measurement uncertainty, but also attempt to include the uncertainty that results from geometry 
imperfections and facility operations. These additional considerations will likely rely on both simulation and 
measured data. A more in-depth analysis of the turbulence models for this massively shock separated 
flow will also be performed. In particular, a detailed comparison of the turbulence energy budgets may 
shed some light on why the shock-separated corner and side wall portions of the flow displayed unchecked 
turbulence production, while other shock-separated regions did not. Finally, the models will be subjected to 
different facility stagnation conditions and back pressures for validation over more than just the target 80% 
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normal shock pressure condition, and to flesh out any Reynolds number effects. 
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